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Abstract 

By using the wavelet transformation (WT), we have analyzed the response 
of an ensemble of A'^ (=1, 10, 100 and 500) Hodgkin-Huxley (HH) neurons to 
transient M-pulse spike trains (M = 1 — 3) with independent Gaussian noises. 
The cross-correlation between the input and output signals is expressed in 
terms of the WT expansion coefficients. The signal-to-noise ratio (SNR) is 
evaluated by using the denoising method within the WT, by which the noise 
contribution is extracted from output signals. Although the response of a 
single (A^ = 1) neuron to sub-threshold transient signals with noises is quite 
unreliable, the transmission fidelity assessed by the cross-correlation and SNR 
is shown to be much improved by increasing the value of N: a population 
of neurons play an indispensable role in the stochastic resonance (SR) for 
transient spike inputs. It is also shown that in a large-scale ensemble, the 
transmission fidelity for supra-threshold transient spikes is not significantly 
degraded by a weak noise which is responsible to SR for sub-threshold inputs. 
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I. INTRODUCTION 



In recent years much studies have been made for the stochastic resonance (SR), in which 
weak input signals are enhanced by background noises 0. This paradoxical SR phe- 
nomenon was first discovered in the context of climate dynamics, and it is now reported 
in many non-linear systems such as electric circuits, ring lasers, semiconductor devices and 
neurons. 

For single neurons, SR has been studied by using various theoretical models such as 
the integrate-and-fire (IF) model the FitzHough-Nagumo (FN) model 0- [§] and 

the Hodgkin-Huxley (HH) model In these studies, a weak periodic (sinusoidal) signal 
is applied to the neuron, and it has been reported that the peak height of the interspike- 
interval (ISI) distribution [Q- or the signal-to-noise ratio (SNR) of output signals 
shows the maximum when the noise intensity is changed. 

SR in coupled or ensemble neurons has been also investigated by using the IF model |T^- 
FN model |T^- |T^ and HH model |T5[- [|T8[. The transmission fidelity is examined by 
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calculating various quantities: the peak-to- width ratio of output signals [jT2| [|T^ [jT6|, 
the cross-correlation between input and output signals [|10| [|18[, SNR |]T0[ |]13[ ||18||, and the 
mutual information . One or some of these quantities has been shown to take a maximum 
as functions of the noise intensity and the coupling strength. Collins, Chow and Imhoff \T} 
have pointed out that SR of ensemble neurons is improved as the size of an ensemble is 
increased. Some physiological experiments support SR in real, biological systems of crayfish 



ig m, cricket |T1 and rat m [|23|. 

SR studies mentioned above are motivated from the fact that peripheral sensory neurons 
play a role of transducers, receiving analog stimula and emitting spikes. In central neural 
systems, however, cortical neurons play a role of data-processors, receiving and transmitting 
spike trains. The possibility of SR in the spike transmission has been reported |]2^- ||26|| . 
The response to periodic coherent spike-train inputs has been shown to be enhanced by 
an addition of weak spike-train noises whose interspike intervals (ISIs) have the Poison or 
gamma distribution. 

It should be stressed that these theoretical studies on SR in neural systems have been 
performed mostly for stationary analog (or spike-train) signals although they are periodic 
or aperiodic There has been few theoretical studies on SR for non-stationary signals. 
Fakir |p8[ has discussed SR for non- stationary analog inputs with finite duration by calculat- 

have 
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ing the cross-correlation. By applying a single impulse, Pei, Wilkinson and Moss 
demonstrated that the spike-timing precision is improved by noises in an ensemble of 1000 
HH neurons. One of the reason why SR study for stationary signals is dominant, is mainly 
due to the fact that the stationary signals can be easily analyzed by the Fourier transforma- 
tion (FT) with which, for example, the SNR is evaluated from FT spectra of output signals. 
The FT requires that a signal to be examined is stationary, not giving the time evolution of 
the frequency pattern. Actual biological signals are, however, not necessarily stationary. It 
has been reported that spike signals in cortical neurons are generally not stationary, rather 
they are transient signals or bursts whereas periodic spikes are found in systems such 
as auditory systems of owl and the electrosensory system of electric fish [pl |. 

The limitation of the FT analysis can be partly resolved by using the short-time Fourier 
transformation (STFT). Assuming that the signal is quasi-stationary in the narrow time 
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period, the FT is applied with time-evolving narrow windows. Then STFT yields the time 
evolution of the frequency spectrum. The STFT, however, has a critical limitation violating 
the uncertainty principle, which asserts that if the window is too narrow, the frequency 
resolution will be poor whereas if the window is too wide, the time resolution will be less 
precise. This limitation becomes serious for signals with much transient components, like 
spike signals. 

The disadvantage of the STFT is overcome in the wavelet transformation (WT) In 
contrast to the FT, the WT offers the two-dimensional expansion for a time-dependent signal 
with the scale and translation parameters which are interpreted physically as the inverse of 
frequency and time, respectively. As a basis of the WT, we employ the mother wavelet which 
is localized in both frequency and time domains. The WT expansion is carried out in terms 
of a family of wavelets which is made by dilation and translation of the mother wavelet. 
The time evolution of frequency pattern can be followed with an optimal time-frequency 
resolution. 

The WT appears to be an ideal tool for analyzing signals of a non-stationary nature. 
In recent years the WT has been applied to an analysis of biological signals |3^, such as 
electoencephalographic (EEG) waves ||3^- ||4^, and spikes |^5|. EEG is a reflection of 



the activity of ensembles of neurons producing oscillations. By using the WT, we can obtain 
the time-dependent decomposition of EEG signals to 6 (0.3-3.5 Hz), 9 (3.5-7.5 Hz), a (7.5- 
12.5 Hz), p 



'12.5-30.0 Hz) and 7 (30-70 Hz) components It has been shown that 

the WT is a powerful tool to the spike sorting in which coherent signals of a single target 
neuron are extracted from mixture of response signals ^l|- [^. Quite recently Hasegawa 



| 45| has made an analysis of transient spike-train signals of a HH neuron with the use of 
WT, calculating the energy distribution and Shanon entropy. 

It is interesting to analyze the response of ensemble neurons to transient spike inputs 
in noisy environment by using the WT. There are several sources of noises: (i) cells in 
sensory neurons are exposed to noises arising from the outer world, (ii) ion channels of the 
membrane of neurons are known to be stochastic 
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m 



the synaptic transmission yields 
noises originating from random fluctuations of the synaptic vesicle release rate [0, and 
(iv) synaptic inputs include leaked currents from neighboring neurons |^8|. Most of existing 
studies on SR adopt the Gaussian noises, taking account of the items (i)-(iii). Simulating 
the noise of the item (iv), Refs. |2^- [^] include spike-train noises whose ISIs have the 
Poisson or gamma distribution. In this study we take into account Gaussian noises; SR for 
spike-train inputs with added spike-train noises will be discussed in a separate paper ^ . 

We assume an ensemble of Ho dgkin- Huxley (HH) neurons to receive transient spike trains 
with independent Gaussian noises. The HH neurons model is adopted because it is known 

The signal transmission is assessed 
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to be the most realistic among theoretical models 
by the cross-correlation between input and output signals and SNR, which are expressed 
in terms of WT expansion coefficients. In calculating the SNR, we adopt the denoising 
technique within the WT method ||5l|- |5^, by which the noise contribution is extracted 
from output signals. 

Our paper is organized as follows: In Sec. HA, an adopted model for an ensemble of 
iV-unit HH neurons is described, and in Sec. IIB the WT is briefly discussed. We present 
the calculated results for M = 1 pulse train in Sec. IIIA, and the results for spike trains 
with M = 2 and 3 are discussed in Sec. IIB. The final Sec. IV is devoted to conclusion and 
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discussion. 



II. CALCULATION METHODS 

A. Ensemble Neuron Model 

We assume a network consisting of A^-unit HH neurons which receive the same spike 
trains but independent Gaussian white noises through excitatory synapses. Spikes emitted 
by the ensemble neurons are collected by a summing neuron. A similar model was previously 
adopted by several authors studying SR for analog signals [|1^] [|T5| [0. An input signal 
in this paper is a transient spike train consisting of M impulses (M = 1 — 3). Dynamics 
of the membrane potential Vi of the HH neuron i is described by the non-linear differential 
equations given by 

C dVi{t)/dt = -/r + /f + /f, (for l<t<N) (1) 

where = 1 /iF/cm^ is the capacity of the membrane. The first term 1^°'^ of Eq.(l) denotes 
the ion current given by 

ir = 9n.m%{V, - ^Na) + QKutiV, - Vk) + - VI), (2) 

where the maximum values of conductivities of Na and K channels and leakage are (^Na = 
120 mS/cm^, g^, = 36 mS/cm^ and gi^ = 0.3 mS/cm^, respectively; the respective reversal 
potentials are = 50 mV, Vk = —77 mV and Vi^ = —54.5 mV. Dynamics of the gating 
variables of Na and K channels, mi, hi and rij, are described by the ordinary differential 



equations, whose details have been given elsewhere pOl p6| 



The second term Jf^ in Eq.(l) denotes the postsynaptic current given by 
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ir= Y.9s{Va-V,)a{t-t,^), (3) 



m=l 



which is induced by an input spike with the magnitude Va given by 



M 



U,{t) = Va E ^(^-^i-), (4) 
m=l 

with the alpha function a{t): 

a{t) = {t/n) e-*/- e(t). (5) 

In Eqs.(3)-(5) is the m-th firing time of the input, the Heaviside function is defined by 
6(t) = 1 for X > and for x < 0, and gs, and stand for the conductance, reversal 
potential and time constant, respectively, of the synapse. 

The third term If in Eq.(l) denotes the Gaussian noise given by 

< >= 0, (6) 



4 



< Jr(t) I^{t') >= 2D 5u Sit - t'), (7) 

where the overhne X and the bracket < X > denote the temporal and spatial averages, 
respectively, and D the intensity of white noises. 

The output spike of the neuron i in an ensemble is given by 

Uoiit) = VaY. ^(t-toin), (8) 

n 

in a similar form as an input spike [Eq.(4)], where tom is the n-th firing time when Vi{t) 
crosses = mV from below. 

We should remark that our model given by Eqs.(l)-(7) does not include couplings among 
ensemble neurons. This is in contrast with some works on ensemble neurons [|l3l |17 



I8| where introduced couplings among neurons play an important role in SR besides noises, 
related discussion being given in Sec. IV. 

Differential equations given by Eqs.(l)-(7) are solved by the forth-order Runge-Kutta 
method by the integration time step of 0.01 ms with double precision. Some results are 
examined by using the exponential method. The initial conditions for the variables are 
given by Vi{t)= -65 mV, mi(t)=0.0526, /ii(t)=0.600, ni{t)=0.313 at t = 0, which are the 
rest-state solution of a single HH neuron. Hereafter time, voltage, conductance, current, 
and D are expressed in units of ms, mV, mS/cm^, /iA/cm^, and /iA^/cm^, respectively. We 
have adopted parameters of K = 30, Vc = —50, and Tg = r„ = 2. Adopted values of gs, D, 
M and will be described shortly. 



B. Wavelet Analysis 

There are two types of WTs: one is the continuous wavelet transformation (CWT) 
and the other the discrete wavelet transformation (DWT). In the former the parameters 
denoting the scale and trasnlation are continuous variables while in the latter they are 
discrete variables. 

The CWT for a given regular function f{t) is defined by 

c(a, b) = Jdt Cfe(t) fit) ^Pabit), fit) >, (9) 
with a family of wavelets ipabit) generated by 

i^abit)=\a\-'/'iji^-^), (10) 
a 

where ipit) is the mother wavelet, the star denotes the complex conjugate, and a and b express 
the scale change and translation, respectively, and they physically stand for the inverse of 
the frequency and the time. Then the CWT transforms the time-dependent function fit) 
into the frequency- and time-dependent function c(a, 6). The mother wavelet is a smooth 
function with good localization in both frequency and time spaces. A wavelet family given by 
Eq.(lO) plays a role of elementary function, representing the function fit) as a superposition 
of wavelets ipabit)- 

The inverse of the wavelet transformation may be given by 
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f{t) = C^'J ^/d6c(a,6)^.,(t), (11) 

when the mother wavelet satisfies the following two conditions: 
(i) the admissibility condition given by 

< Cy, < oo, (12) 

with 

/oo 
duj \^{uj) / \uj \, (13) 



where ^(a;) is the Fourier transformation of ip{t), and 
(ii) the zero mean of the mother wavelet: 



dt il^it) = ^(0) = 0. (14) 



On the contrary, the DWT is defined for discrete values of a = 2^ and h = 2% (j, k: 
integers) as 

Cjk = c(2^ 2^k) =< i/jjkit), fit) >, (15) 

with 

^jkit)=2~^^^^lj{2-H-k). (16) 
The ortho-normal condition for the wavelet functions is given by 

< '^jk{t), l/^j'k'it) >= Sjf Skk': (17) 

which leads to the inverse DWT: 

/W = E E CjkMt)- (18) 

j k 

In the multiresolution analysis (MRA) of the DWT, we introduce a scaling function (j){t), 
which satisfies the recurrent relation with 2K masking coefficients, hk, given by 

2K-1 

(t>{t) = ^2 J2 hkcP{2t-k), (19) 

k=0 

with the normalization condition for (j){t) given by 

J dt 0(i) = 1. (20) 
A family of wavelet functions is generated by 

2K-1 

^{t) = V2 ^ (-l)^/i2i^_i_fc <P{2t - k). (21) 

fe=0 



6 



The scaling and wavelet functions satisfy the ortho-normal relations: 



< (Pit), (Pit -m) > =5^0, (22) 
<^it),^ljit-m) > =6mo, (23) 

< 0(t),V^(t-m) > = 0. (24) 

A set of masking coefficients hj is chosen so as to satisfy the conditions shown above. 

The simplest wavelet function for K = 1 is the Harr wavelet for which we get ho = hi = 
l/v/2, and 

ipnit) = 1, for < t < 1/2 
= -1, for 1/2 < t < 1 

= 0, otherwise (25) 

In the more sophisticated wavelets like the Daubechies wavelet, an additional condition given 

by 

J dtt^ i^it) = 0, for £ = 0,1, 2, 3 L-1 (26) 

is imposed for the smoothness of the wavelet function. Furthermore, in the Coiflet wavelet, 
for example, a similar smoothing condition is imposed also for the scaling function as 

J dtt^ (Pit) = 0, for £ = 1,2, 3 L'-l (27) 

Once WT coefficients are obtained, we can calculate various quantities such as auto- 
and cross-correlations and SNR, as will be discussed shortly. In principle the expansion 
coefficients Cjk in DWT may be calculated by using Eqs.(15) and (16) for a given function 
fit) and an adopted mother wavelet ipit)- This integration is, however, inconvenient, and in 
an actual fast wavelet transformation, the expansion coefficients are obtained by a matrix 
multiplication with the use of the iterative formulae given by the masking coefficients and 
expansion coefficients of the neighboring levels of indices, j and k p2 . 

One of the advantages of the WT over FT is that we can choose a proper mother 
wavelet among many mother wavelets, depending on a signal to be analyzed. Among many 
candidates of mother wavelets, we have adopted the Coiflet with level 3, compromising the 
accuracy and the computational effort. The WT has been performed by using the MATLAB 
wavelet tool box. 



III. CALCULATED RESULTS 

A. Input Pulses with M = 1 

Firstly we discuss the case in which ensemble HH neurons receive a common single 
(M = 1) impulse. When input synaptic strength is small: Qs < Qth, no neurons flre in the 
noise-free case, while it is sufficiently large igg > Qth) neurons flre, where gth = 0.091 is the 
threshold value. For a while, we will discuss the sub-threshold case of Qs = 0.06 < gth with 
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= 500. The M = 1 input pulse is applied at t = 100 ms, as shown in Fig. 1(a). Figure 
1(b) shows the time dependence of the postsynaptic current Jj = Jf^ + If of the neuron 
i = 1 with added noises of D = 0.5. Because neurons receive noises for < t < 100 ms, 



the states of neurons when they receive input pulse are randomized [|16[. Figure 1(c) shows 
the time dependence of the membrane potential Vi of the neuron 1, which fires with a delay 
of about 6 ms. This delay is much larger than the conventional value of 2-3 ms for the 
supra-threshold inputs because the integration of marginal inputs at synapses needs a fairly 



a long period before firing [^. Firings in ensemble 500 neurons for D = 0.5 are depicted 
by raster in Fig. 2(a). We note that neurons fire not only by input pulses plus noises but 
also spuriously by noises only. When the noise intensity is increased to D = 1.0, spurious 
firings are much increased, as Fig. 2(b) shows. 

We assume that information is carried by firing times of spikes but not by details of their 
shapes. In order to study how information is transmitted through ensemble HH neurons with 
the use of the DWT, we divide the time scale by the width of time bin of Tb as t = = £ Tb 
{£: integer), and define the input and output signals within the each time bin by 

M 

Wi{t) = J2 e(|t-ti„|-Tb/2), (28) 

m=l 



N 

Woit) = (1/iV) 0(1 t - I -Tb/2). (29) 

i=l n 

In Eqs. (28) and (29) Q{t) stands for the Heaviside function, Wi{t) the external input signal 
(without noises), Wo(t) the output signal averaged over the ensemble neurons, the m-th 
firing time of inputs, and tom the n-th firing time of outputs of the neuron i [Eq.(8)]. The 
time bin is chosen as Tb =5 ms in our simulations. A single simulation has been performed 
for 320 (=2^ Tb) ms. Figures 1(d) and 1(e) show the time dependence of Wi(t) and Wo(t), 
respectively, for the case of D = 0.5. The magnitude of Wo{t) is much smaller than that 
of Wi(t) because only a few neurons fire among 500 neurons. The peak position of Wo(t) is 
slightly shifted compared with that of Wi{t) because of a significant delay of neuron firings 
as mentioned above. 

Wavelet Transformation 

Now we apply the DWT to input and output signals . By substituting f(t) = Wi(t) or 
Wo(t) in Eq.(15), we get their WT coefficients given by 

cxjk = Jdt r.kit) Wx{t), (A = i, o) (30) 

where ipjk{t) is a family of wavelets generated from the mother Coiflet wavelet [see Eq.(16)]. 
The uppermost frame of Fig. 3(a) expresses the input signal Wi{t). Note that the lower and 
upper horizontal scales express b and b Tt (in units of ms), respectively. Figure 3(b) shows the 
calculated WT coefficients of Wi{t) which are plotted as a function of b{j) = {k — 0.5) 2^ for 
various j values after convention. The WT decomposition of the signal: f{t) = Z]j=i fjit), 
is plotted in Fig. 3(a). The WT coefficients of j = 1 and 2 have large values near b = 
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20 {b Tb = 100 ms) where Wi{t) has a peak. Contributions from j = 1 and j = 2 are 
predominant in Wi{t). It is noted that the WT coefficient for j = 4: has a significant value at 
6 ~ 56 far away from b = 20. The WT decomposition and WT coefficients for output signal 
Wo{t) for D = 0.5 are shown in Fig. 4(a) and 4(b), respectively. The dominant contribution 
arises from j = 1 in Wo{t). Similar plots of the WT decomposition and the WT coefficients 
of the output signal Woit) for D = 1.0 are presented in Figs. 5(a) and 5(b), respectively. As 
the noise intensity is increased, fine structures in the WT coefficients appear, in particular 
for small j. 

Auto- and Cross- correlations 

The auto-correlation functions for input and output signals are defined by 

Taa = J dt Wxity Wx{t), (31) 

= M-'Y.Y. ^Ijk CA.fc, (A= i, o) (32) 

j k 

where the ortho- normal relations of the wavelets given by Eqs.(22)-(24) are employed. Sim- 
ilarly the cross-correlation between input and output signals is defined by 

TUf3) = J dt W,{ty Wo{t + /3 Tb), (33) 
= ^"' E E c*fcCo,fc(/?), (34) 

j k 

where c^k and CojkiP) are the expansion coefficients of Wi{t) and Wo(t + l3Ti,), respectively. 
The maxima in the cross-correlation and the normalized one are given by 

T = maXf3[TioiP)], (35) 



= maxp[ J^°^^^ ]. (36) 

It is noted that for the suprashreshold inputs in the noise-free case, we get Ta = Too = Tio = 1 
and then F = 7 = 1. 

Figure 6(a) shows the D dependence of F, Fqo and 7 for (yf = 0.06 and = 500. They 
are zero at D = because of the adopted sub-threshold parameters. As increasing D from 
zero, F and Fqo are gradually increased. Because of the factor of I/a/Fqo in Eq.(36), the 
magnitude of 7 is larger than those of F and Fqo- We note that 7 is enhanced by weak noises 
and it is decreased at larger noises, which is a typical SR behavior. 

Signal-to-Noise Ratio 



We will evaluate the SNR by employing the denoising method [Q. The key point in 
the denoising is how to choose which wavelet coefficients are correlated with the signal and 
which ones with noises. The simple denoising method is to neglect some DWT expansion 
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coefficients when reproducing the signal by the inverse wavelet transformation. We get the 
denoising signal by the inverse WT [Eq.(18)]: 



wtit) = T.T.^t;kMtl (37) 

j k 

with the desnoising WT coefficients cjl to be chosen properly as will be discussed below. 
The simplest denoising method, for example, is to assume that WT components for a < Oc 
in the (a, h) plane arise from noises to set the denoising WT coefficients as 

= Cjfc, for i > jc (a > Oc) 

= 0, otherwise (38) 



where jc (= logg etc) is the critical j value |^ 



In this study we adopt a more sophisticated method, assuming that the components for 
6 < 6l or 6 > 6u at a < ac in the (a,b) plane are noises to set the denoising WT coefficients 

as 

c'jk = Cjk, for j > jc or ki,<k <ku 

= 0, otherwise (39) 

where jc (= log2 Cc) denotes the critical j value, and /cl (= &l S^-') and k\j (= bu 2~^) are 
the lower and upper critical k values, respectively. We get the inverse, denoising signal by 
using Eq.(37) with the denoising WT coefficients determined by Eq.(39). 

From the above consideration, we may define the signal and noise components by 

As = Jdt W^^'ity W^'^it), (40) 

= EEIco^l^ (41) 



An = jdt [WM Wo{t) - <"(t)* <"(t)], (42) 

=i:i:i\<^o,k?-\c%\\ (43) 

3 k 

The SNR is defined by 

SNR = lQ\og,^{AJAn) (dB). (44) 

In the present case we can fortunately obtain the WT coefficients for ideal case of noise- 
free and supra-threshold inputs. We then properly determine the denoising parameters of 
jci hi and hjj. From the observation of the WT coefficients for the ideal case, which is not 
shown here but is not dissimilar to those shown in Figs. 3(a) and 4(a), we assume that the 
upper and lower bounds, may be chosen as 

&L = toi/Tb - 5h, 6u = W/Tb + 6h, (45) 
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where toi (^om) are the first (M-th) firing times of output signals in the ideal case of noise- 
free and supra-threshold inputs, and Sb denotes the marginal distance from the b values 
expected to be responsible to the signal transmission. 

The denoising WT coefficients c^^ of Wo for D = 0.5 and D = 1.0 are denoted by dashed 
bars in Fig. 4(b) and 5(b), respectively. Note that the positions of dashed bars are slightly 
shifted leftward to avoid the superposition with the solid bars representing the original WT 
coefficients. In the case oi D = 0.5, the denoising WT coefficients are almost the same as 
the original ones. The inverse WT signal of PF^^", which is given by Eq.(37) and which is 
shown by the dashed curve in Fig. 4(a), is almost the same as the original signal shown by 
the solid curve. On the contrary, in the case of D = 1.0, the denoising WT coefficients of 
Wo do not include the terms far away from b = 20 for j < jc in the original WT coefficients. 
The inverse WT signal of Wo^, shown by the dashed curve in Fig. 5(a), has less amount of 
ripples than the original signal shown by the solid curve. 

Figure 7(a) and 7(b) show the calculated SNR of Wo as a function of 6b for various jc 
values with D = 0.5 and D = 1.0, respectively. We note that the value of SNR is rather 
insensitive to a choice of the parameters of 6b and jc- Then we have decided to adopt 6b = 5 
and jc = 3 for our simulations. 

Figure 6(b) shows the D dependence of SNR calculated for g = 0.06 and = 500. We 
note that SNR shows a typical SR behavior: a rapid rise to a clear peak and a slow decrease 
for larger value of D. 



So far our simulation has been made for the fixed parameters of gs = 0.06 and N = 500, 
which will be changed in the followings. Figure 8(a) and 8(b) show the D dependence of 
the cross-correlation 7 and SNR, respectively, calculated for = 1, 10 100 and 500. The 
results for = 1 and 10 are averaged values of ten and hundred trials, respectively, their 
root-mean-square (RMS) values being shown by error bars in Fig. 8(a) and 8(b). The results 
of N = 500 show the typical SR behavior as mentioned before [Figs. 6(a) and 6(b)]. On the 
contrary, SR effect for a single {N = 1) is marginally realized in SNR but not in 7. Large 
error bars for the results of iV = 1 implies that the reliability of information transmission is 
very low in the sub-threshold condition This is clearly seen in Fig. 9(a) and 9(b) where 
7 and SNR for = 1 are plotted against the trial number of 100 runs. The signal can be 
transmitted only when the signal fortunately coincides with noises and the signal plus noise 
crosses the threshold level. Then only 13 among 100 trials are succeeded in the transmission 
of M = 1 inputs through a neuron. We note in Figs. 8(a) and 8(b) that as the size of the 
network is much increased, the peak values of 7 and SNR are much enhanced and the SR 
behavior becomes more evident. Figures 8(a) and 8(b) demonstrate that the ensemble of 
neurons play an indispensable role for information transmission of transient spike signals. 
This is consistent with the results of Collins et al. [0] and Pei et al. [|l^, who have pointed 



out the improvement of the information transmission by increasing the size of the network. 

Next we change the value of gs, the strength of input synapses. The gg dependence of 
7 and SNR for A^ = 500 neurons is shown in Figs. 10(a) and 10(b), respectively, where 
values oi D = 0.0, 0.5 and 1.0 are employed. Note that in the noise-free case {D = 0), 
we get 7 = 1 and SNR=oo for gs > gth but 7 = and SNR=— cxd for gs < gth, as shown 
by dotted curves in Figs. 10(a) and 10(b). It is shown that moderate sub-threshold noises 
considerably improve the transition fidelity. We also note that the presence of such noises 
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does not significantly degrade the transmission fidelity for supra-thresliold cases in ensemble 
neurons. For a comparison, the Qs dependence of 7 and SNR for single (A^ = 1) neurons is 
shown in Figs. 11(a) and 11(b), respectively. As having been shown in Fig. 8(a), SR for 
iV = 1 is less significant compared to that for N = 500. 



B. Input Pulses with M = 2 and 3 

Now we discuss the cases of M = 2 and 3. Input pulses are applied at t = 100 and 125 
ms for the M = 2 case. The input ISI is assume to be = 25 ms because spikes with this 



value of ISI are reported to be ubiquitous in cortical brains ||29|. Firings of 500 neurons 
for the noise intensity of D = 1.0 are shown by raster in Figs. 12(a), which shows that 
firings occur mainly at t ~ 103 and 128 ms with a delay of about 3 ms. The output signal 
Woit) averaged over = 500 neurons is depicted by the solid curve in the uppermost frame 
of Fig. 13(a), which has main two peaks. The WT decomposition and WT coefficients of 
Wo{t) are plotted in Figs. 13(a) and 13(b), respectively. When we compare the results for 
M = 2 shown in Figs. 13(a) and 13(b) with those for M = 1 shown in Fig. 5(a) and 
5(b), we note that components for j = 1 and 2 for M = 2 are much increased because of 
the presence of the second peak for M = 2 inputs while the contributions from j > 3 are 
changed httle. The denoising has been made by the procedure given by Eqs. (37), (39) 
and (45). The denoising WT coefficients are plotted by dashed bars in Fig. 13(b) and the 
denoising component signals are shown by dashed curves in Fig. 13(a). 

We apply the M = 3 pulse at t = 100, 125 and 150 ms, the input ISI being again Ti = 25 
ms. Raster in Fig. 12(b) shows firing of 500 neurons for D = 1.0. Firings mainly occur at 
t ~ 103, 128 and 153 ms with a delay of 3 ms. The solid curve in the uppermost frame of 
Fig. 14(a) shows the averaged output of Wo(t) for D = 1.0. The sold bars in Fig. 14(b) 
show its WT coefficients and the solid curves in Fig. 14(a) express its WT decomposition. 
Again components of j = 1 and 2 are increased because of the second and third peaks in 
Wo{t). The denoising WT coefficients are shown by dashed bars in Fig. 14(b) and denoising 
signals by dashed curves in Fig. 14(a). 

The calculated D dependence of the cross-correlation 7 [SNR] for M = 1, 2 and 3 is 
plotted in Fig. 15(a) [15(b)]. Both 7 and SNR show typical SR behavior irrespective of the 
value of M, although a slight difference exists between the M dependence of 7 and SNR: 
for larger M, the former is larger but the latter is smaller at the moderate noise intensity 
of D < 1.0. When similar simulations are performed for different ISI values of Ti = 15 and 
35 ms, we obtain results which are almost the same as that for Tj = 25 (not shown). This 
is because the output spikes for inputs with M = 2 and 3 are superposition of an output 
spike for a M = 1 input when the ISI is larger than the refractory period of neurons. 



IV. CONCLUSION AND DISCUSSION 

It has been controversial how neurons communicate information by action potentials or 
spikes [|5^- |6^. The one issue is whether information is encoded in the average firing rate 
of neurons {rate code) or in the precise firing times {temporal code). Since Andrian [K 



first noted the relationship between neural firing rate and stimulus intensity, the rate-code 
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model has been supported in many experiments of motor and sensory neurons. In the last 
several years, however, experimental evidences have been accumulated, indicating a use of 
the temporal code in neural systems [^. Human visual systems, for example, have 
shown to classify patterns within 250 ms despite the fact that at least ten synaptic stages 
are involved from retina to the temporal brain The transmission times between two 
successive stages of synaptic transmission are suggested to be no more than 10 ms on the 
average. This period is too short to allow rates to be determined accurately. 

Although much of debates on the nature of the neural code has focused on the rate versus 
temporal codes, there is the other important issue to consider: information is encoded in the 
activity of single (or very few) neurons or that of a large number of neurons {population or 
ensemble code). The population rate code model assumes that information is coded in the 
relative firing rates of ensemble neurons, and has been adopted in the most of the theoretical 
analysis. On the contrary, in the population temporal code model, it is assumed that relative 
timings between spikes in ensemble neurons may be used as an encoding mechanism for 
perceptional processing |7^- [|72[. A number of experimental data supporting this code 



have been reported in recent years |7^- ||76[. For example, data have demonstrated that 
temporally coordinated spikes can systematically signal sensory object feature, even in the 
absence of changes in firing rate of the spikes [ [F^ . 

Our simulations based on the temporal-code model has shown that a population of 
neurons plays a very important role for the transmission of sub-threshold transient spike 
signals [Figs. 8(a) and 18(b)]. In particular for single neurons the transmission is quite 
unreliable and the appreciable SR effect is not realized. When the size of ensemble neurons 
is increased, the transmission fidelity is much improved in a fairly wide-range of parameters 
including both the sub- and supra-threshold cases (Fig. 10). We note in Figs. 8(a) and 
8(b) that 7 (or SNR) for = 100 is different from and larger than that for A^ = 1 with 100 
trials. This seems strange because a simulation for an ensemble of A^ neurons is expected 
to be equivalent to simulations for a single neuron with A^ trials, if there is no couplings 
among neurons as in our model. This is, however, not true, and it will be understood as 
follows. We consider a quantity of X{N, A^) which is 7 (or SNR) averaged over A^^ trials 
for an ensemble of A^ neurons. We implicitly express X{N, Nj.) as 



X{N, Nr) = < F(< >) > (46) 

1 ^ 1 ^ U)-, 



E^^E^r), (47) 



with 



W-, = w 



W = E 0(1^-^1^ I -V2), (48) 



where -C ■ ^ and < ■ > stand for averages over trials and an ensemble neurons, respectively, 
defined by Eq.(47), t^'^^ is the n-th firing time of the neuron i in the /x-th trial, wl^\t) is 
its output signal of the neuron i within a time bin of Tb [cf. Eq.(29)], and F{y{t)) is a 
functional of a given function of y{t) relevant to a calculation of 7 (or SNR). Figures 8(a) 
and 8(b) show that the relation: X(100, 1) > X(l, 100), namely 

F{< ^ >) X F{w[^^) >, (49) 
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holds for our 7 and SNR. Note that if F(-) is hnear, we get X(100, 1) = X(l, 100). This 
imphes that the inequahty given by Eq. (49) is expected to arise from a non-/mear character 
of F{-). This reminds us of the algebraic inequality: f{< x >) > < f{x) > valid for a 
convex function f{x), where the bracket < ■ > stands for an average over a distribution 
of a variable x. It should be again noted that there is no couplings among our neurons 
in the adopted model. Then the enhancement of SNR with increasing N is only due to a 
population of neurons. This is quite different from the result of some papers |]T3 |T7 



p!8| in which the transmission fidelity is enhanced not only by noises but also by introduced 
couplings among neurons in an ensemble. 

In our simulations reported in Sec. Ill, independent noises are applied to ensemble 
neurons. If instead we try to apply the same or completely correlated noise to them, it 
is equivalent to applying noises to a single neuron, and then appreciable SR effect is not 
reahzed as discussed above. Then SR for transient spikes requires independent noises to 
be applied to a large-scale ensemble of neurons. This is consistent with the result of Liu, 



Hu and Wang [jT8[ who discussed the effect of correlated noises on SR for stationary analog 
inputs. 

Although spike trains with small values of M = 1 — 3 have been examined in Sec. Ill, we 
can make an analysis of spikes with larger M or bursts, by using our method. In order to 
demonstrate its feasibility, we have made simulations for transient spikes with larger M. The 
upper curve of Fig. 16(a) expresses input Wi{t) with a M = 7 spike train whose firing times 
are ti^= 100, 115, 130, 145, 160, 180, and 200 ms, and ISIs are T; =15 and 20 ms. Firings 
of 100 neurons in an ensemble is depicted by raster in Fig. 16(b), where the parameters 
of Qs = 0.06, D = 1.0 and N = 100 are adopted. The lower curve in Fig. 16(a) shows 
its output Wo(t) averaged over the ensemble. We apply the WT to Wo(t) to get its WT 
coefficients and its WT decomposition, the latter being shown in Fig. 16(b). The j = 1 and 
j = 2 components are dominant. After the denoising, we get 7 = 0.523 and S'A^i?=18.6 dB, 
which are comparable to those for D = 1.0 with M = 1 — 3 shown in Figs. 15(a) and 15(b). 
In our denosiging method given by Eqs.(39) and (45), we extract noises outside the b region 
relevant to a cluster of spikes, but do not take account of noises between pulses. When a 
number of pulses M and/or the ISI become larger, a contribution from noises between 
pulses become considerable, and it is necessary to modify the denoising method such as to 
extract noises between pulses, for example, as given by 

Cjfc = Cjk, for j > jc or /cl™ <k< k\jm (m = 1 - M) 

= 0. otherwise (50) 

In Eq. (50) kim and are m- and j-dependent lower and upper limits given by 

A;Lm = 2-^' (tom/Tb - 5h), fcun, = 2"^' {t^^/T^ + 5h) , (51) 

where tom is the m-th firing time for the noise-free and supra-threshold input and 5h the 
margin of h. 

An ensemble of neurons adopted in our model can be regarded as the front layer (referred 



to as the layer I here) of a synfire chain [77]: output spikes of the layer I are feed- forwarded 
to neurons on the next layer (referred to as the layer II) and spikes propagate through the 
synfire chain. The postsynaptic current of a neuron I on the layer II is given by 
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N 



^/N) E E (Va - Vs) ait - re, - t,,m) + (52) 



where wa and tu are the synaptic couphng and delay time, respectively, for spikes to propa- 
gate from neuron i on the layer I to neuron C. on layer II, toim the m-th firing time of neuron 
i on the layer I, and noises. Transmission of spikes through the synfire chain depends 
on the couplings (excitatory or inhibitory) wu, the delay time r^j and noises I^. There are 
some discussions on the stability of spike transmission in a synfire chain. Quite recently 
Diesmann has shown by simulations using IF models that the spike propagation with a 
precise timing is possible in fairly large-scale synfire chains with moderate noises. 

By augmenting our neuron model including the coupling term given by Eq.(52) and 
tentatively assuming wm = w = l.b and = 0, we have calculated the transmission fidelity 
by applying our WT analysis to output signals on the layer II as well as those on the layer 
I. Calculated cross-correlation and SNR are shown in Fig. 17(a) and 17(b), respectively. In 
Fig. 17(a), 7l (L=I, II) denotes the cross-correlation between the input and output signals 
on the layer L. Similarly, SNRl in Fig. 17(b) expresses the SNR of output signal on the 
layer L. Note that 7/ and SNRi are nothing but 7 and SNR having been shown in Figs. 
6(a) and 6(b). We note that the transmission fidelity on the layer II is better that that of the 
layer I because 7/ < 7// and SNRj < SNRjj at almost D values. We have chosen w = 1.5 
such that neurons on the layer II can fire when more than 6% of neurons on the layer I fire. 
When we adopt smaller (larger) value of w, both 7// and SNRjj abruptly increase at larger 
(smaller) D value. However, the general behavior of the D dependence of 7// and SNRji is 
not changed. This improvement of the transmission fidelity in the layer II than in the front 
layer I is expected to be partly due to the threshold-crossing behavior of HH neurons. It 
would be interesting to investigate the transmission of signals in a synfire chain by including 
SR of its front layer. 

Our paper entirely relies on numerical simulations. We are currently trying to work on 
theoretical description of the result reported in this paper. Conventional approaches having 
been employed for a study of SR such as the rate-equation and linear-response theories 
|T|]- 0], do not work on our case. Mato |^ adopted Gammaitoni's approach []7S[ for an 



analysis of his SR result with the continuous spike-train signals. It seems, however, not to 
be transposed directly to our case of transient spike-train signals even if our HH model is 
replaced by simpler IF model or threshold-crossing model. Its analytical study is left as our 
future problem. 

To summarize, the response of an ensemble of neurons to transient spike trains has been 
discussed with the use of the WT. Although the transmission of sub-threshold transient 
inputs is not reliable in a single neuron, it is much improved in a large-scale ensemble of 
neurons with a weak noise (SR effect). It is also shown that this noise does not significantly 
degrade the transmission fidelity for supra-threshold inputs in the large-scale ensemble. A 
population of neurons plays an essential role in SR for transient inputs. 
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FIGURES 



FIG. 1. The time dependence of (a) the M = 1 input pulse (Ui), (b) the postsynaptic current 
(Ii = /f^ + If), (c) the membrane potential (Vi), and (d) input (Wi) and (e) output signals (Wq) 
within the time bin of T^, = 5 ms for D = 0.5, gg = 0.06 and N = 500. 

FIG. 2. Rasters showing firings in ensemble neurons for (a) D = 0.5 and (b) D = 1.0 with 
M = 1, gs = 0.06 and N = 500. 

FIG. 3. (a) The input signal Wi(t) (uppermost frame) for M = 1, and its WT decomposition: 
/ = Z)j=i fj- (b) Its WT expansion coefficients Cjk- Curves of fj for j > 2 and of / in (a), and 
those of Cjk for j > 2 in (b) are successively shifted upward by 1.0. Upper horizontal scale expresses 
b Tb in units of ms. 

FIG. 4. (a) The output signal Wo{t) (uppermost frame) for D = 0.5 and M = 1, and its 
WT decomposition: / = J2j=i /?• (b) 1^ WT expansion coefficients Cjk- Solid and dashed curves 
denote the original and denoising results, respectively. Same as in Fig. 3. 

FIG. 5. (a) The output signal Wo{t) (uppermost frame) for D = 1.0 and M = 1, and its 
WT decomposition: / = Ej=i fj- (b) Its WT expansion coefficients Cjk- Solid and dashed curves 
denote the original and denoising results, respectively. Same as in Fig. 3. 

FIG. 6. The D dependence of 7 (circles), Too (triangles) and F (squares) for M = 1, gg = 0.06 
and = 500. 

FIG. 7. The 5b dependence of SNR of output signals for (a) D = 0.5 and (b) D = 1.0 with 
jc = 2 (circles), 3 (triangles) and 4 (squares) for M = 1, gg = 0.06 and N = 500, the denoising 
condition being given by Eqs.(39) and (45). 

FIG. 8. (a) The cross-correlation 7 and (b) SNR as a function of D for M = 1 and gg = 0.06 

with = 1, 10, 100 and 500. Results for iV = 1 (A^ = 10) arc averaged values of 100 (10) trials and 
their RMS arc shown by error bars. RcsTilts of (a) and (b) for N = 10, 100 and 500 are successively 
shifted upward by 1.0 and 20, respectively. 

FIG. 9. (a) The cross-correlation 7 and (b) SNR of 100 trials for single (A'' = 1) neurons with 
gg = 0.06 and M = 1. 

FIG. 10. The gg dependence of the cross-correlation 7 and SNR for N = 500 and M = 1, 
Dot, dashed and solid curves denoting the results of D = 0.0, 0.5 and 1.0, respectively. The arrow 
expresses the threshold gg value. 
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FIG. 11. The Qs dependence of the cross-correlation 7 and SNR for A'^ = 1 and M = 1: Dot, 
dashed and sohd curves denote the results of -D = 0.0, 0.5 and 1.0, respectively, and error bars 
express RMS values for 100 trials. 

FIG. 12. Raster showing firings in ensemble neurons for (a) M = 2 and (b) M = 3 with 
D = 1.0, gs = 0.06 and N = 500. 

FIG. 13. (a) The output signal Wo{t) (uppermost frame) for D = 1.0 and M = 2, and its 
WT decomposition: / = Ej=i fj- (b) Its WT expansion coefficients Cjk- Solid and dashed curves 
denote the original and denoising results, respectively. Same as in Fig. 3. 

FIG. 14. (a) The output signal Wo{t) (uppermost frame) for D = 1.0 and M = 3, and its 
WT decomposition: / = J2j=i fj- (b) Its WT expansion coefficients Cj^- Solid and dashed curves 
denote the original and denoising results, respectively. Same as in Fig. 3. 

FIG. 15. (a) The cross-correlation 7 and (b) SNR as a function of D for M = 1 (circles), 
M = 2 (squares) and M = 3 (triangles). 

FIG. 16. (a) The input spike train Wi{t) (upper curve) consisting of M = 10 pulses, and its 
output signal Wo{t) (lower curve) for gg = 0.06, D = 1.0 and N = 100. (b) Raster showing firings 
of 100 neurons, (c) The WT decomposition of Wo{t), results for j =2, 3, 4 and 5 being successively 
shifted upward by 0.2. 

FIG. 17. (a) The D dependence of the cross-correlation 7/ (filled circles) and 7// (open circles), 
and (b) that of SNRj (filled circles) and SNRu (open circles) for M = 1; note that SNRjj = 00 
for 0.3 < L> < 0.6 (see text). 
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